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Energy and X^ax reconstruction of hadron-initiated showers in surface arrays 
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Abstract: The current methods to determine the primary energy in surface arrays are different when deaUng with hadron 
or photon initiated showers. In this work, we adapt a method previously developed for photon-initiated showers to hadron 
primaries. We determine the Monte Carlo parametrizations that relate the surface energy estimator and the maximum of 
shower development for both, proton and Iron primaries. Using for each primary their own set of calibration curves, 
which is of course impossible in practice, we show that the energy could be inferred with a negligible bias and 12% 
resolution. However, we show that a mixed calibration could also be performed, including both type of primaries, such 
that the bias still remains low and the achieved resolution is around 15%. In addition, the method allows the simultaneous 
determination of Xmax in pure surface arrays with resolution better than 7%. 
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1 Introduction 

In large ultra-high energy cosmic ray surface arrays the 
technique traditionally used to determine the primary en- 
ergy consists on the inference of the lateral distribution 
of particles of the extensive air shower (EAS) it pro- 
duces into the atmosphere. Scintillators (e.g., Volcano 
Ranch, AGASA, KASCADE) and water Cherenkov detec- 
tors (e.g., Haverah Park) have been mainly used for this 
purpose. The surface array allows the discrete sampling of 
the shower front at ground level and then, the lateral dis- 
tribution is fitted assuming a certain functional form (LDF, 
lateral distribution function). Later, the signal inferred at 
an optimum distance is used as the energy estimator, which 
is related to the primary energy using Monte Carlo (MC) 
parametrizations. This optimum distance is traditionally 
fix for each experiment since it is assumed to be only de- 
pendent on the array spacing. However, recent studies sug- 
gest the convenience of calculating the optimum distance 
for each individual shower taking into account primary en- 
ergy and direction [2] . A special case from the experimen- 
tal point of view is the Pierre Auger Observatory which pi- 
oneers the simultaneous use of water Cherenkov detectors 
and fluorescence telescopes. For these hybrid events, sys- 
tematic errors in their energy estimate are greatly reduced 
[3]. 

These methods assume that primary is a proton and it is 
considered adequate for heavier primaries since the esti- 
mated energy for nuclei depends weakly on their mass 
number. On the other hand, the difference is signifi- 



cant when dealing with photon-initiated showers since the 
muonic component is much lower, shower development is 
affected by the geomagnetic field and the LPM effect de- 
lays in average the first interaction. Thus, each experiment 
has followed a different method when searching for photon 
primaries. Haverah Park and Auger use the muon density 
at ground [4] and MC parametrizations [5] to infer the pri- 
mary energy respectively, while in [6, 7], AGASA data is 
directly compared to photon MC simulations. 

The method used by Auger [5] was first proposed in [1]. 
The original idea was to adapt the energy reconstruction 
method to the late development of photon showers and the 
key point is to rely explicitly on the development stage of 
the shower. The method was originally applied to photon 
showers in an Auger-like array. An empirical parametriza- 
tion between 6'(1000)/£^ and AX = Xground/ cos9 - 
Xmax was found, where 5'(1000) is the inferred signal at 
1000 meters from the shower axis, E is the primary en- 
ergy, Xmax is the maximum of shower development and 
is the zenith angle of the shower. The parametrization is 
almost independent of the primary energy due to the well- 
known universality of the electromagnetic component of 
EAS [8, 9, 10] and the small muon fraction in the photon- 
initiated showers. 

In this work, we show how to modify this method to be ap- 
plicable to hadron-initiated showers where the muon com- 
ponent is significant, specially, in case of water Cherenkov 
arrays which enhanced their contribution to the total mea- 
sured signal. In this way it is possible to determine simul- 
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taneously the energy and X^ax of the shower from pure 
surface data. Alternatively, several surface parameters have 
been used to infer indirectly Xmax^ such as the rise time of 
the signals in the detectors and the azimuthal features of 
the time distributions [11]. 

2 Simulations 

The simulation of the atmospheric showers is performed 
with the AIRES Monte Carlo program (version 2.8.4a) 
[12] using QGSJET-II-03 as the hadronic interaction model 
(HIM). The input primary energy goes from log{E/eV) = 
19.0 to 19.6 in 0.1 steps. Around 300 events have been 
simulated per each energy and for both, proton and iron 
primaries. The zenith angle has been selected following a 
sine-cosine distribution from 30 to 60 degrees, while the 
azimuth angle is randomly distributed. 

As it will be explained in Section 3, only the reconstructed 
5(1000) and zenith angle of the shower will be needed in 
the method proposed here. We simulate the reconstructed 
zenith angle by fluctuating the real one with a Gaussian 
whose standard deviation is 1^, a typical value of its uncer- 
tainty in surface arrays [13, 14, 15, 16]. 

A realistic S'(IOOO) of the event could be obtained from 



5'(1000)(i9) 



X [l^Cx- Dx^] 



(1) 



where x = cos^{9) — cos^ (38^). A, B, C and D are con- 
stants obtained in [17] for QGSJet-II-03 and for iron and 
proton primaries. However, event by event fluctuations and 
reconstruction uncertainties are not taken into account if 
5(1000) is directly estimated from Eq. (1). Thus, we cal- 
culate the reconstructed 5(1000) from our own simulation 
of the detector response and reconstruction process, as it is 
explained next. 

The simulation of the tank response and the fit of the LDF 
are performed using our own simulation program previ- 
ously tested in [2]. An infinite array whose unitary cell 
consists on a triangle with detectors separated 1.5 km is 
considered. The real core is randomly located inside an el- 
ementary cell and the reconstructed one is determined fluc- 
tuating the real one with a Gaussian function whose stan- 
dard deviation depends on the primary energy and compo- 
sition (more details en [2]). 

The signal collected at each station for a given shower with 
a certain energy and zenith angle is set assuming a true 
lateral distribution function like [17] 



5(r) = 5(1000) x(^)- 



r + 700 



(2) 

^ 1000 + 700 ' ' ^ ^ 

where the distance to shower axis r is in meters and 
/3((9, 5(1000)) is given by 

a -\- b{secO — 1) ^ if secO < 1.55 
f3{0, 5(1000)) = <( a + blsecO - 1) + f{secO - 1)^, 

if secO > 1.55 

(3) 



where a = 2.26 + 0.195/o^(e), b = -0.98, c = 0.37 - 
OMsecO + O.SOsec'^O, d = 1.27 - 0.27sec0 + O.OSsec'^O, 
e = c(5(1000))^ and / = -0.9. 5(1000) in Eqs. (2) 
and (3) is obtained from Eq. (1). Later, the signal assigned 
to each station is fluctuated using a Poissonian distribution 
whose mean is given by the true LDF. We impose typical 
values for trigger condition and saturation. Saturated de- 
tectors are excluded from the LDF fit. 

Next, the lateral distribution of particles is fitted using a 
functional form given by 



log 5(r) = ai + a2 



r + 700 
1000 + 700 



(4) 

where the slope of the LDF and the normalization constant 
are determined in each fit considering the core position as 
fixed in the reconstructed one. Finally, the reconstructed 
5(1000) is determined as the interpolated value from the 
fit at 1000 meters from the shower axis. 



3 Energy and X^ax reconstruction 

The method is based on MC parametrizations and, essen- 
tially, the energy of the primary and X^nax are iteratively 
determined. As mentioned, the only ingredients needed are 
the reconstructed zenith angle of the incoming shower, 6, 
and the interpolated 5(1000) value from the LDF fit. 

The iterative process should be started with an initial rough 
estimation of the energy. We arbitrarily use 10, 20 and 50 
EeV. The results do not depend on this choice. Next, Xmax 
is estimated using its average dependence on energy. 



Xrr 



(7o + X logio{Eprim) g/cm^, (5) 



where Eprim is in EeV. Next, using 5(1000) and 6, primary 
energy could be obtained from 



5(1000)/^pWm=2^0X 



1 



AX-lOO 

Pi 



l + ( 



P2 



VEM/EeV. 

(6) 

Then, Eprim obtained from Eq. (6), is used in Eq. (5) to 
get a new Xmax following an iterative process until con- 
vergence is achieved (we set the convergence in energy at 
10~^). The convergence is always fast, only 3-6 steps are 
needed. 

This method works properly in case of photon-initiated 
showers since S {!{){){))/ Eprim vs. AX shows an univer- 
sal profile independently of primary energy due to the low 
fraction of muons in photon showers as explained before. 
However, hadron primaries produce EAS with a significant 
muon component and, in addition, the muon fraction de- 
pends on the primary energy. As consequence, if a global 
calibration curve is used, the inferred energy will show an 
energy dependent bias as shown in Fig. 1 (black points). 

In order to account for the muonic component properly, dif- 
ferent fits could be performed for each primary energy (Fig. 
2). Then, the reconstructed energy obtained with the global 
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Figure 1 : Energy error vs. primary energy for proton show- 
ers using their own caHbration. The points and the error 
bars are the mean and sigma of the Gaussian fit to the en- 
ergy error distribution respectively. Similar results for iron 
primaries. 
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Figure 2: 6'(1000)/£^ vs AX for proton showers. The fits 
corresponds to Eq. (6). Similar curves are obtained for iron 
showers. 

fit could be used to select a new calibration curve from Fig. 
2, repeating the process until the reconstructed energy is 
not modified. Only one or two iterations are needed. Thus, 
the energy bias is corrected (Fig. 1, red points) while res- 
olution remains below 12% considering the whole simula- 
tions set. 

On the other hand, Xmax is also properly reconstructed. In 
fact, the bias is negligible and resolution is around 7% for 
proton and 4% for iron primaries respectively. 

Obviously, it is impossible in practice to use different cal- 
ibration curves for each primary. However, we could also 
calculate a mixed calibration fitting both type of primaries 
simultaneously (Fig. 3). The energy iX^nax) error is shown 
in Fig. 4 (Fig. 5). It remains below 10% (7%) for each 
primary and it is almost negligible if both are considered 
together. Resolution is better than 15% (7%). 



Figure 3: As Fig. 2 but mixing iron and proton primaries. 
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Figure 4: Energy error vs. primary energy for proton, iron 
and both primaries using the mixed calibration. The points 
and the error bars are the mean and sigma of the Gaussian 
fit to the energy error distribution respectively. 
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Figure 5: As Fig. 4 but for X„ 



Finally, we have determined the energy and Xmax biases in 
case of a mixed composition sample. We have selected 100 
samples with 100 events each. Proton and Iron primaries 
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are randomly included in the sample such the proton frac- 
tion varies from to 1 in 0.1 steps. The energy {Xmax) bias 
as a function of the proton fraction is shown in Fig. 6 (Fig. 
7). Using the mixed calibration, the energy {Xmax) bias 
varies from —7% (—5%) to +6% (+5%) as composition 
changes from pure proton to pure Iron. 
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Figure 6: Energy bias as a function of the proton fraction 
using the three sets of calibration curves analyzed in this 
work. The points and the error bars are the median and the 
confidence levels at 68 and 95% respectively. 



in the inferred energy varies from —7% to +6% as compo- 
sition changes from pure proton to pure Iron, while resolu- 
tion is better than 15%. 

In addition, the method allows the indirect determination 
of the maximum of shower development, Xmax, from pure 
surface data. The resolution achieved is around 7% and the 
bias goes from —5% for pure proton to +5% in case of pure 
Iron. Therefore, Xmax could be estimated from surface 
data reliably, whose statistics are around 10 times larger 
than that from fluorescence telescopes. However, it is im- 
portant to note that MC parametrizations could be affected 
by the fact that the simulations do not reproduce properly 
experimental data [18]. 
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Figure 7: As Fig. 6 but for X„ 



4 Conclusions 

An iterative method, previously developed to infer the pri- 
mary energy of photon-induced showers in pure surface ar- 
rays, has been modified to be applicable to hadron-initiated 
showers. The inferred energy bias is negligible and the 
resolution is around 12% if each primary could be recon- 
structed with its own calibration. In a more realistic ap- 
proach, the set of calibration curves have been also ob- 
tained mixing proton and Iron primaries and then, the bias 



[I] P Billoir et al. arXiv:astro-ph/0701583, 2007. 

[2] G. Ros et al.NIM A. 608, Issue 3, p. 454-463, 2009. 
[3] The Pierre Auger Collaboration. Proc. 31th ICRC. 

Lodz, Poland, arXiv:0906.2189, 2009 
[4] M. Ave et al. Phys.Rev. D65, 063007, 2002. 
[5] Pierre Auger Collaboration. Astroparticle Physics, 29, 

p. 243-256, 2008. 
[6] M. Risse et al. PRL 95, 171102, 2005. 
[7] G.I. Rubtsov et al. Phys.Rev.D 73:063009, 2006. 
[8] S. Lafebre et al. arXiv:0902.0548, 2009. 
[9] W. D. Apel et al. Astropart. Phys. 29:412-419, 2008. 
[10] F. Schmidt et al. arXiv: 0706. 1990, 2007. 

[II] The Pierre Auger Coflaboration. Proc. 31th ICRC. 
Lodz, Poland, 2009. arXiv:0906.2319, 2009. 

[12] S. Sciutto. http://www.fisica.unlp.edu.ar/auger/aires. 
[13] The Pierre Auger Collaboration. Proceedings of the 

29th ICRC, Pune, V. 7, p. 17, 2005. 
[14] AGASA Collaboration. Astropart.Phys. 10, 303-311, 

1999. 

[15] HiRes Coflaboration. Nucl. Phys. Proc. Suppl. 165, 

239-245, 2007. 
[16] M. Ave et al. Astropart. Phys. 19, 47-60, 2003. 
[17] loana C. Maris. PhD Thesis, http://digbib.ubka.uni- 

karlsruhe.de/volltexte/documents/697773 
[18] The Pierre Auger Coflaboration. arXiv:0906.2319, 

2009. 



